Nonlinear and chaotic ice ages: data vs speculations 
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It is shown that, the wavelet regression detrended fluctuations of the reconstructed temperature 
for the past 400,000 years (Antarctic ice cores data) are completely dominated by one-third subhar- 
monic resonance, presumably related to Earth precession effect on the energy that the intertropical 
regions receive from the Sun. Effects of Galactic turbulence on the temperature fluctuations are 
also discussed. Direct evidence of chaotic response of the atmospheric CO2 dynamics to obliquity 
periodic forcing has been found in a reconstruction of atmospheric CO2 data (deep sea proxies), for 
the past 650,000 years. 

PACS numbers: 92.70.Gt, 92.70.Qr, 92.10.am, 92.10.Hm 
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Forcing of nonlinear systems does not always have the 
results that one might expect. Subharmonic and chaotic 
resonances are prominent examples of such response. Re- 
cent paleoclimate reconstructions provide indications of 
nonlinear properties of Earth climate at the late Pleis- 
tocene [H,[l| (the period from 0.8 Myr to present). Long 
term decrease in atmospheric CO2, which could result 
in a change in the internal response of the global car- 
bon cycle to the obliquity forcing, has been mentioned as 
one of the principal reasons for this phenomenon (see, for 
instance, Q-Q). At present time one can recognize at 
least three problems of the nonlinear paleoclimate, which 
we will address in present paper using recent data and 
speculations. 

A. Reconstructed air temperature on millennial time 
scales are known to be strongly fluctuating. See, for 
instance figures 1 and 2. While the nature of the 
trend is widely discussed (in relation to the glaciation 
cycles) the nature of these strong fluctuations is still 
quite obscure. The problem has also a technical as- 
pect: detrending is a difficult task for such strong fluctu- 
ations. In order to solve this problem a wavelet regres- 
sion detrending method was used in present investiga- 
tion. Then a spectral analysis of the detrended data re- 
veals rather surprising nature of the strong temperature 
fluctuations. Namely, the detrended fluctuations of the 
reconstructed temperature are completely dominated by 
so-called one-third subharmonic resonance, presumably 
related to Earth precession effect on equatorial insola- 
tion. 

B. Influence of Galactic turbulent processes on the 
Earth climate can be very significant for time-scales less 
than 2.5 kyr. 

C. Nonlinearity of the Earth climate can also results 
in a chaotic response to an external periodic forcing. Al- 
ready pioneering studies of the effect of external peri- 
odic forcing on the first Lorenz model of the chaotic cli- 
mate (convection) revealed very interestingproperties of 
chaotic response (see, for instance, [(| , [?] , 8] ) . The cli- 
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FIG. 1: The reconstructed air temperature data (dashed line 
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(see also Ref. [121). The solid curve (trend) corresponds to a 
wavelet (symmlet) regression of the data. 



mate, where the chaotic behavior was discovered for the 
first time, is still one of the most challenging areas for 
the chaotic response theory. One should discriminate be- 
tween chaotic weather (time scales up to several weeks) 
and a more long-term climate variation. The weather 
chaotic behavior usually can be directly related to chaotic 
convection (as it was done for the first time by Lorenz), 
while appearance of the chaotic properties for the long- 
term climate events is a non-trivial and challenging phe- 
nomenon. It was suggested that such properties can play 
a significant role for glaciation cycles at multi-millennium 
time scales |3|.j^|.(l0]. 

Cyclic forcing, due to astronomical modulations of the 
solar input, rightfully plays a central role in the long- 
term climate models. Paradoxically, it is a very non- 
trivial task to find imprints of this forcing in the long- 
term climate data. It will be shown in present paper that 
just unusual properties of nonlinear and chaotic response 
are the main source of this problem. 
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FIG. 2: The wavelet regression detrended fluctuations from 
the data shown in Fig. 1. 
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FIG. 3: Spectrum of the wavelet regression detrended fluctu- 
ations shown in Fig. 2. 
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Figure 1 shows reconstructed air temperature data 
(dashed line) for the period 0-340 kyr as presented at the 
site (Antarctic ice cores data, see also Ref. 12]). The 
solid curve (trend) in the figure corresponds to a wavelet 
(symmlet) regression of the data (cf Ref. [Hj]). Fig- 
ure 2 shows corresponding detrended fluctuations, which 
produce a statistically stationary set of data. Most of 
the regression methods are linear in responses. At the 
nonlinear nonparametric wavelet regression one chooses 
a relatively small number of wavelet coefficients to rep- 
resent the underlying regression function. A threshold 
method is used to keep or kill the wavelet coefficients. 
In this case, in particular, the Universal (VisuShrink) 
thresholding rule with a soft thresholding function was 
used. At the wavelet regression the demands to smooth- 
ness of the function being estimated are relaxed consid- 
erably in comparison to the traditional methods. Figure 
3 shows a spectrum of the wavelet regression detrended 
data calculated using the maximum entropy method (be- 
cause it provides an optimal spectral resolution even for 
small data sets). One can see in this figure a small peak 
corresponding to period ~ 5kyr and a huge well defined 
peak corresponding to period ^15kyr. We also obtained 
analogous results (approximately 10% larger) from the 
"Vostok" ice core data for period 0-420kyr (for the data 
description see Refs. [3,[l|). 

In order to understand underlying physics of the very 
characteristic picture shown in the Fig. 3 (cf. Ref. [lit ]) 
let us imagine a forced excitable system with a large 
amount of loosely coupled degrees of freedom schemat- 
ically represented by Duffing oscillators (which has be- 
come a classic model for analysis of nonlinear phenomena 
and can exhibit both deterministic and chaotic behavior 
[l7|-(20| depending on the parameters range) with a wide 
range of the natural frequencies loq : 



where x denotes the temporal derivative of x, ft is the 
strength of nonlinearity, and F and lo are characteristic 
of a driving force. It is known (see for instance Ref. [20j p 
that when lo ~ 3lo and ft <C 1 the equation (1) has a 
resonant solution 
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where the amplitude a and the phase tp are certain con- 
stants. This is so-called one-third subharmonic reso- 
nance with the driving frequency lo corresponding ap- 
proximately to 5kyr period (the huge peak in Fig. 3 
correspond to the first term in the right-hand side of 
the Eq. (2)). For the considered system of the oscilla- 
tors an effect of synchronization can take place and, as 
a consequence of this synchronization, the characteristic 
peaks in the spectra of partial oscillations coincide plj ]. 
It can be useful to note, for the climate modeling, that 
the odd-term subharmonic resonance is a consequence of 
the reflection symmetry of the natural nonlinear oscilla- 
tors (invariance to the transformation x — > —x). 

Origin of the periodic energy input with the period ~ 
5kyr can be related to dynamics of the energy that the 
intertropical regions receive from the Su n ( equatorial in- 
solation). Indeed, it is found in Ref. [22| that a clear 
and significant 5kyr period is present in this dynamic 
over last 1 Myr. The amplitude of the 5kyr cycle in the 
insolation decreases rapidly when getting away from the 
equator. Using the fact that double insolation maximum 
and minimum arise in the tropical regions in the course 
of one year, the authors of the Ref. [22j speculated that 
this period in seasonal amplitude of equatorial insolation 
is determined by fourth harmonic of the Earth precession 
cycle. It should be noted, that the idea of a significant 
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FIG. 4: A small-time-scales part of the autocorrelation func- 
tion defect of the wavelet regression detrended fluctuations 
from the data shown in Fig. 2. The straight line indicates 
the Kolmogorov's '2/3' power law for the structure function 
in ln-ln scales. 



role of tropics in generating long-term climatic variations 
is rather a new one (see Ref. |22j | for relevant references). 
In Ref. [23J , for instance, the authors speculated that the 
high frequency climate variability (in the millennial time 
scales) could be related to high sensitivity of the trop- 
ics to summer time insolation. Then, the oceanic ad- 
vective transport could transmit an amplified response 
of tropical precipitation and temperature to high lati- 
tudes. Physical mechanism of this amplification is still 
not clear and the above discussed one-third subharmonic 
resonance can be a plausible possibility (in this respect 
it is significant, that we used the Antarctic data). 



GALACTIC TURBULENCE AND THE 
TEMPERATURE FLUCTUATIONS 



Since the high frequency part of the spectrum is cor- 
rupted by strong fluctuations (the Nyquist frequency 
equals 0.5 ](500y) -1 ]), it is interesting to look at cor- 
responding autocorrelation function C(t) in order to un- 
derstand what happens on the millennial time scales. 
Figure 4 shows a relatively small-times part of the cor- 
relation function defect. The ln-ln scales have been used 
in this figure in order to show a power law (the straight 
line) for structure function: {(v(t + r) — v(t)) 2 ) : 



1 - C(t) oc ((v(t + r) - v(t)) 2 ) oc r 2 / 3 



(3) 



Although, the scaling interval is short, the value of 
the exponent is rather intriguing. This exponent is well 
known in the theory of fluid (plasma) turbulence and cor- 
responds to so-called Kolmogorov's cascade process. This 
process is very universal for turbulent fluids and plasmas 

m 



271 ] . For turbulent processes on Earth and in Helio- 



This power law: '2/3', for structure function (by virtue 
of the Taylor hypothesis transforming the time scaling 
into the space one (2^].[25|) is known for fully developed 
turbulence as Kolmogorov's power law. 



sphere the Kolmogorov-like scaling with such large time 
scales certainly cannot exist. Therefore, one should think 
about a Galactic origin of Kolmogorov turbulence with 
such large time scales (let us recall that diameter of the 
Galaxy is approximately 100,000 light years). This is not 
surprising if we recall possible role of the galactic cosmic 
rays for Earth climate (see, for instance, [10], [28]- 30]). 
In this respect, it should be also noted that the '2/3' 
scaling law can reflect not the velocity field of the galac- 
tic interstellar media but the turbulent electron density 
field (presumably produced by supernova) according to 
the Obukhov-Corrsin turbulent mixing process (see, for 
instance, Ref. [HI], [III and references therein). 



ATMOSPHERIC C0 2 DYNAMICS AT 
MULTI-MILLENNIUM TIME SCALES 



The angle between Earth's rotational axis and the 
normal to the plane of its orbit (known as obliquity) 
varies periodically between 22.1 degrees and 24.5 de- 
grees on about 41, 000- year cycle. Such multi- millennium 
timescale changes in orientation change the amount of 
solar radiation reaching the Earth in different latitudes. 
In high latitudes the annual mean insolation (incident so- 
lar radiation) decreases with obliquity, while it increases 
in lower latitudes. Obliquity forcing effect is maximum 
at the poles and comparatively small in the tropics. Mi- 
lankovic theory suggests that lower obliquity, leading to 
reduction in summer insolation and in the mean insola- 
tion in high latitudes, favors gradual accumulation of ice 
and snow leading to formation of an ice sheet. The obliq- 
uity forcing on Earth climate is considered as the primary 
driving force for the cycles of glaciation (see for a recent 
review [l[). Observations show that glacial changes from 
-1.5 to -2.5 Myr (early Pleistocene) were dominated by 
41 kyr cycle 0,(13], j34|, whereas the period from 0.8 Myr 
to present (late Pleistocene) is characterized by approx- 
imately 100 kyr glacial cycles [35j],[36|]. While the 41 kyr 
cycle of early Pleistocene glaciation is readily related to 
the 41 kyr period of Earth's obliquity variations the 100 
kyr period of the glacial cycles in late Pleistocene still 
presents a serious problem. Influence of the obliquity 
variations on global climate started amplifying around 
2.5 Myr, and became nonlinear at the late Pleistocene. 
Long term decrease in atmospheric CO2, which could 
result in a change in the internal response of the global 
carbon cycle to the obliquity forcing, has been mentioned 
as one of the principal reasons for this phenomenon (see, 
for instance, Therefore, investigation of the his- 
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FIG. 5: A reconstruction of atmospheric CO2 based on deep- 
sea proxies, for the past 650kyr. The data were taken from 

M. 
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FIG. 6: Spectrum of atmospheric CO2 fluctuations for the 
data shown in Fig. 5 



toric variability in atmospheric CO2 can be crucial for 
understanding the global climate changes at millennial 
timescales. 

Figure 5 shows a reconstruction of atmospheric CO2 
based on deep-sea proxies, for the past 650kyr (the data 
taken from |37[). Resolution of the data set is 2kyr. 
Fluctuations with time-scales less than 2kyr could be 
rather large (statistically up to 308ppm [37(), but they 
are smoothed by the resolution. Figure 6 shows a power 
spectrum of the data set calculated using the maximum 
entropy method. The spectrum exhibits a peak indi- 
cating a periodic component (the arrow in the Fig. 6 
indicates a lOOkyr period) and a broad-band part with 



A semi-logarithmic plot was used in Fig. 6 in order to 
show the exponential decay more clearly (at this plot the 
exponential decay corresponds to a straight line). Both 
stochastic and deterministic processes can result in the 
broad-band part of the spectrum, but the decay in the 
spectral power is different for the two cases. The ex- 
ponential decay indicates that the broad-band spectrum 
for these data arises from a deterministic rather than a 
stochastic process. For a wide class of deterministic sys- 
tems a broad-band spectrum with exponential decay is a 
generic feature of their chaotic solutions Refs. [38j-[43j|. 
Let us recall that the discussed above Duffing oscillators 
can exhibit both deterministic and chaotic behavior (l7j - 
[2£ ' 



38J depending on the parameters range. 



Nature of the exponential decay of the power spectra 
of the chaotic systems is still an unsolved mathematical 
problem. A progress in solution of this problem has been 
achieved by the use of the analytical continuation of the 
equations in the complex domain (see, for instance, [43| ) . 
In this approach the exponential decay of chaotic spec- 
trum is related to a singularity in the plane of complex 
time, which lies nearest to the real axis. Distance be- 
tween this singularity and the real axis determines the 
rate of the exponential decay. For many interesting cases 
chaotic solutions are analytic in a finite strip around the 
real time axis. This takes place, for instance for attrac- 
tors bounded in the real domain (the Lorentz attractor, 
for instance). In this case the radius of convergence of 
the Taylor series is also bounded (uniformly) at any real 
time. If parameters of the dynamical system fluctuate pe- 
riodically around their mean values with period T e , then 
the restriction of the Taylor series convergence (at certain 
conditions) is determined by the period of the parametric 
modulation, and the width of the analytic strip around 
real time axis equals T e /2n [n|. Let us consider, for sim- 
plicity, solution u(t) with simple poles only, and to define 
the Fourier transform as follows 

u(w) = (2^)~ 1 / 2 [ T °' 2 dt e~ lut u{t) (5) 

J-Te/2 

Then using the theorem of residues 

u(oj) = z(27r) 1 / 2 Rj exp(iwxj — \tJyj\) (6) 
3 

where Rj are the poles residue and Xj + iyj are their lo- 
cation in the relevant half plane, one obtains asymptotic 
behavior of the spectrum E(lj) — |u(w)| 2 at large u> 



E(u) ~ exp(-2|cjy min |) 



(7) 



where y m i n is the imaginary part of the location of the 
pole which lies nearest to the real axis. Therefore, expo- 
nential decay rate of the broad-band part of the system 
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spectrum, Eq. (4), equals the period of the parametric 
forcing. 

The chaotic spectrum provides two different character- 
istic time-scales for the system: a period corresponding to 
fundamental frequency of the system, Tt un , and a period 
corresponding to the exponential decay rate, T e = l// e 
(cf. Eq. (4)). The fundamental period Tf un can be esti- 
mated using position of the low- frequency peak, while the 
exponential decay rate period T e = l// e can be estimated 
using the slope of the straight line of the broad-band part 
of the spectrum in the semi-logarithmic representation 
(Fig. 6). From Fig. 6 we obtain Tf un ~ 95 ± 8 kyr (the 
peak is quite broad due to small data set) and T e ~ 41 ±1 
kyr (the estimated errors are statistical ones). 

Thus, the obliquity period of 41 kyr is still a dominat- 
ing factor in the chaotic CO2 fluctuations, although it is 
hidden for linear interpretation of the power spectrum. 
In the nonlinear interpretation the additional period 
Tf un ~ 100 kyr may correspond to the fundamental 
frequency of the underlying nonlinear dynamical system 
and it determines the apparent 100 kyr 'periodicity' 
of the glaciation cycles for the last 650 kyr (cf Refs. 
0:05 44 1 and references therein). 



The data were provided by World Data Center for 
Paleoclimatology, Boulder and NOAA Paleoclimatology 
Program. I also acknowledge that a software provided by 
K. Yoshioka was used at the computations. 
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